function [beta_J_hat, R2, beta_J_hat_tod, R2_tod,tod] = Jump_Beta(data, factors, L_k, B_kl, un)

% This code allows nan for missing data, but require qn to be divisor of n
% for speed and better finite sample performance.

% L_k: the input should be a list, with number of jump factors for each
% factor; sum of L_k = H
% B_kl: the input should be a H*2 matrix. Specifying the range of the H
% jump factors.
% data: stock return table {intra_time, permno, ret, prc, date, symbool}
% factors: [date, seconds, factors 1-6, risk free rate]

    
    M         =      1;
    K         =      size(L_k,1);
    H         =      sum(L_k);
       
    stock_dates  = year(data{:,5})*10000 + month(data{:,5})*100 + day(data{:,5});
    dates        = unique(stock_dates);
    
    factor_stock = [factors(:,1), factors(:,3:2+K), factors(:,end), NaN(size(factors,1),1)];
    
    for j = 1:size(dates)
        factor_stock(factors(:,1)==dates(j),end) = data{stock_dates==dates(j),3};
    end
    
    nDays     =      size(unique(factor_stock(:,1)),1);
    n         =      size(factor_stock,1);
    nobs      =      n/nDays;
    Delta_n   =      1/252/(nobs-1);
    Tn        =      n*Delta_n;
    vn        =      1000*log(1/(Tn*Delta_n));
    
    dFt       =      factor_stock(:,2:1+K)';
    dFt_bar   =      zeros(H, n);
    
    dPt       =      factor_stock(:,end)' - factor_stock(:,end-1)'; % Ret - RF
    dPt_zero  =      dPt;
    dPt_zero(isnan(dPt)) = 0;

       
    % Jump regression
    % Use this matrix to filter factor jump to make sure we are considering
    % time period corresponding to each individual stock's live time. 
    dFt_filter_P = zeros(H, n, M);
    for i = 1:M
        dFt_filter_P(:,:,i) = repmat(~isnan(dPt(i,:)), H, 1);
    end
    
    % Step 1 Stack all the factor jump together into dFt_bar(H*n)
    [dFt_cont_ix_tod,dFt_cont_ix,tod] = fn_cont_ix_tod(un, 0.47, dFt, Delta_n, nDays, nobs, []);
    dFt_jp_ix = ones(K,n) - dFt_cont_ix;
    dFt_jp_ix_tod = ones(K,n) - dFt_cont_ix_tod;
    
    id = zeros(H,n);
    id_tod = zeros(H,n);
        
    for i = 1:K
        for j = 1:L_k(i)
            m = sum(L_k(1:i-1));
            id(m+j, :) = dFt_jp_ix(i,:) & dFt(i,:)>min(B_kl(m+j,:)) & dFt(i,:)<max(B_kl(m+j,:)); 
            id_tod(m+j, :) = dFt_jp_ix_tod(i,:) & dFt(i,:)>min(B_kl(m+j,:)) & dFt(i,:)<max(B_kl(m+j,:)); 
        end
    end
    
    for i = 1:K
        for j = 1:L_k(i)
            m = sum(L_k(1:i-1));
            dFt_bar(m+j,:) = dFt(i,:);
        end
    end
    
    dFt_bar = dFt_bar.*id;
    dFt_bar_tod = dFt_bar.*id_tod;

    beta_J_hat = zeros(M,H);
    beta_J_hat_tod = zeros(M,H);
    for i = 1:M
        if min(eig((dFt_bar.*dFt_filter_P(:,:,i))*(dFt_bar.*dFt_filter_P(:,:,i))')) > 1/vn
            beta_J_hat(i, :) = dPt_zero(i,:)*(dFt_bar') *...
                            inv((dFt_bar.*dFt_filter_P(:,:,i))*(dFt_bar.*dFt_filter_P(:,:,i))');
        else
            beta_J_hat(i, :) = 0;
        end
    end

    for i = 1:M
        if min(eig((dFt_bar_tod.*dFt_filter_P(:,:,i))*(dFt_bar_tod.*dFt_filter_P(:,:,i))')) > 1/vn
            beta_J_hat_tod(i, :) = dPt_zero(i,:)*(dFt_bar_tod') *...
                            inv((dFt_bar_tod.*dFt_filter_P(:,:,i))*(dFt_bar_tod.*dFt_filter_P(:,:,i))');
        else
            beta_J_hat_tod(i, :) = 0;
        end
    end
    
    R2         = (beta_J_hat * (dFt_bar.*dFt_filter_P(:,:,1))) * (beta_J_hat * (dFt_bar.*dFt_filter_P(:,:,1)))';
    R2         = R2/((dPt_zero.*(sum(dFt_bar~=0)~=0))*(dPt_zero.*(sum(dFt_bar~=0)~=0))');

    R2_tod         = (beta_J_hat * (dFt_bar_tod.*dFt_filter_P(:,:,1))) * (beta_J_hat_tod * (dFt_bar_tod.*dFt_filter_P(:,:,1)))';
    R2_tod         = R2_tod/((dPt_zero.*(sum(dFt_bar_tod~=0)~=0))*(dPt_zero.*(sum(dFt_bar_tod~=0)~=0))');
    
 
end